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Abstract: This paper presents a generalized framework for segmenting 
closed-contour anatomical and pathological features using graph theory and 
dynamic programming (GTDP). More specifically, the GTDP method 
previously developed for quantifying retinal and corneal layer thicknesses is 
extended to segment objects such as cells and cysts. The presented 
technique relies on a transform that maps closed-contour features in the 
Cartesian domain into lines in the quasi-polar domain. The features of 
interest are then segmented as layers via GTDP. Application of this method 
to segment closed-contour features in several ophthalmic image types is 
shown. Quantitative validation experiments for retinal pigmented epithelium 
cell segmentation in confocal fluorescence microscopy images attests to the 
accuracy of the presented technique. 
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1. Introduction 

Fast, accurate, and objective detection and quantification of imaging biomarkers are crucial 
for the study and diagnosis of ophthalmic diseases. Due to the time consuming and subjective 
nature of manual segmentation, considerable work has been done in recent years to automate 
the segmentation of ocular structures. For example, many layer segmentation algorithms have 
been developed to delineate the retinal [1-6], choroid-scleral [7,8], and corneal layers [9,10] 
on optical coherence tomography (OCT) images [11]. 

However, in addition to layered structures, there is also a need to automatically segment 
and quantify closed-contour features in ophthalmic images. Examples of closed-contour 
anatomical and pathological structures include retinal pigmented epithelial (RPE) cells in 
confocal fluorescence microscopy images of flat-mounted retina [12], intra-retinal cysts in 
OCT images of patients with diabetic macular edema (DME) [13,14], and photoreceptors in 
adaptive optics scanning ophthalmoscope (AOSO) images of retina [15-18]. Previous works 
have been presented to segment closed-contour features, including fluid or detachments on 
OCT images of the retina [19-21], drusen in fundus photographs [22,23] or in OCT images 
[24], geographic atrophy in fundus auto-fluorescence images [25], vessel lumens in 
intravascular OCT images [26], corneal cells [27-31], and photoreceptors on AOSO images 
[32-37]. However, these segmentation algorithms were often developed for specific 
ophthalmic applications, rather than as a general framework applicable to both layer and 
closed-contour features. Furthermore, in some cases such as with photoreceptor or corneal cell 
identification, cell counts or densities were achieved but the actual structure size and shape 
were not necessarily determined [27,29,31-35,37]. 

To segment circular biological structures, transformation of images into the polar domain 
has proven to be an effective tool within the medical community [38-41]. This is due to the 
observation that circles in the Cartesian domain are represented as lines in the polar domain. 
Using this transform, the mathematical operations and manipulations required to segment the 
image are greatly simplified since these operations no longer need to be performed radially. 
However, for non-circular, closed-contour structures, the polar transform is no longer 
adequate since oblong or convoluted features deviate from a flat path. 

In this paper, we present a general framework to segment arbitrarily shaped, closed- 
contour features. We previously developed a framework for layer segmentation based on 
graph theory and dynamic programming (GTDP) [6] to segment retinal and corneal layers on 
spectral domain optical coherence tomography (SDOCT) images [5,6,10]. Section 2 of this 
paper extends the application of our layer segmentation technique to include closed-contour 
segmentation. The algorithm is based on the projection of closed-contour features in the 
Cartesian domain into lines in a virtual domain, which we term the quasi-polar domain. In the 
quasi-polar domain, features of interest are then segmented as layers using our GTDP 
technique. In Section 3, we apply the framework to segment RPE cells in confocal 
nucroscopy images of an age-related macular degeneration (AMD) mouse model, and in 
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Section 4.1 we evaluate the algorithm's effectiveness in detail. Lastly, to highlight the 
flexibihty of the presented framework, Section 4.2 illustrates the segmentation of closed- 
contour features for other ophthalmic features, including cysts on SDOCT and cone 
photoreceptor cells on AOSO images. 

2. A generalized closed-contour segmentation framework using GTDP 

This section introduces a generalized framework for closed-contour feature segmentation. The 
method is an extension of our previously presented technique for layer segmentation via 
GTDP [6]. The core steps are outlined in Fig. 1 and discussed in detail in the following 
subsections. Section 3 then provides the specific algorithmic steps required to implement this 
framework for RPE cell segmentation in confocal fluorescence microscopy images. 



Transform image to I \. Segment layered | \. Transform layer back 
quasi-polar domain structure using GTDP ]—/ to Cartesian domain 



Fig. 1 . Schematic for segmenting closed-contour features via GTDP. 



Generate pilot 
estimate of the 
structure 




2.7. Pilot structure estimation 

The first step is to attain pilot estimates of all structures to segment. These estimates should 
contain information about both the location and shape of the structures. Structure localization 
prevents the algorithm from confusing the target structure with similarly- shaped extraneous 
features, while shape estimation results in a more accurate quasi-polar transform, as will be 
discussed in Section 2.2. 

For a grayscale image if g ^^^""^ in the Cartesian domain, the pilot estimates can be 

represented as a binary image Bf , where 

„ fl, if z e pilot estimate 

Bf(z) = \^ , . ' Vz (1) 

[0, otherwise 

and z = (/, j) denotes the pixel at the ith row and jth column of the image. Many basic 

techniques such as thresholding, contours, local minima, or the multiscale generalized 
likelihood ratio test [42] can be used to attain this estimate, and Section 3.2 describes one such 
method for pilot structure estimation. To isolate a single structure, if and Bf can be 
cropped to smaller image blocks (/^ and BJ, both e 9^^^^ where M <F and N <G. 

2.2. Image transformation into the quasi-polar domain 

This section explains the quasi-polar transform in depth. In its simplest form, the quasi-polar 
transform can be described by a polar transform followed by a flattening step. Figure 2 shows 
an illustration, where a circle, an oval, and an arbitrary shape in the Cartesian domain (Figs. 
2(a-c)) are first transformed into the polar domain (Figs. 2(d-f)). Information about the shape 
of the pilot estimate is then used to flatten the structure into a line in the quasi-polar domain 
(Fig. 2(g)). Once transformed into the quasi-polar domain, the structure can be segmented as 
if it were a layer using the GTDP technique. 

It can also be observed from Fig. 2 that a simple polar transform is sufficient to map 
closed-contour shapes (Figs. 2(a-c)) into layered structures (Figs. 2(d-f)). Motivation for a 
further transformation into the quasi-polar domain (Fig. 2(g)) stems from the tendency for 
shorter geometric paths to be found when using the GTDP technique, especially in low signal- 
to-noise ratio (SNR) imaging scenarios. This is due to the observation that a cut 
(segmentation) with fewer traversed nodes often yields a lower total weight. As a result, 
oblong features are disadvantaged since their paths in the polar domain are longer than the 
shortest geometric path (Figs. 2(e, f)). By performing an extra flattening step, the path of the 
oblong structure now reflects the shortest geometric distance. 
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Fig. 2. Schematic of the quasi-polar transform, (a-c) A circle, oval, and arbitrary shape in the 
Cartesian domain, (d-f) Transformation of (a-c) into the polar domain based on the centroid 
(yellow asterisks in a-c). The centroid in (a-c) becomes a Hne (yellow) in the polar domain, (g) 
Flattening of (d-f) into a Hne in the quasi-polar domain. The result is a transformation of any 
closed-contour shape in the Cartesian domain into a Hne in the quasi-polar domain. 

The quasi-polar transform is implemented using the following three steps: 

1. Map pixels from the Cartesian domain into the polar domain based on their distance 
and angle from a reference pixel and axis, respectively. The reference pixel z^^f can 

be any single pixel, where (z^^f ) = 1; however ideally z^^f lies at the center of the 
closed-contour feature to facilitate its flatness in the polar domain. The centroid of 
the pilot estimate is therefore a good choice for z^^f , where 



1 ^ 

ref ^ (Kef ' J ref ) ^ 



(2) 



for the set of ^pixels {z^lf^p where B^(Z}^) = 1. Example binary images (BJ are 
shown in Figs. 2(a-c) with z^^f marked as a yellow asterisk and the reference axis 
defined as the 7-axis. Using the Cartesian binary image B^ and Eq. (3), generate the 
polar-transformed binary image B^ (Figs. 2(d-f)), where B^(r,0) denotes the pixel 
in B^ with a radius r and angle 0 from the reference pixel and axis, respectively. 
Then, assuming a unit step size, B^ e dl^""^ where R = max(r) and 0 = 361 , 

B^(r,0) = B^(iJ), 

0 < r < max ( J( / - /^^P ' + ( 7 - ) ' ) , V/,7 
0° < 6> < 360°, 

1 = r' sin(6>) + i^^f , 7 = r • cos(0) + j^^^ . 



(3) 



2. Find a function r = f{6) that best estimates the boundary between the background and 
the pilot estimate in B^ . This can be achieved by taking the vertical gradient of the 
image, or by using the GTDP layer segmentation technique described in Section 2.3 
with edge weights determined by the vertical gradient of B^ . 
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3. Generate the quasi-polar binary image e ^i^""® (Fig. 2(g)), where for all columns 
{^pk iLi ' column ft^^ in is determined by 




where C e 



,0x0 



and 




(4) 



0, otherwise 



Use the exact transformation mapping B^ B^ (Steps 2 and Eq. (4)) to then 



Since Step 3 flattens the pilot estimate instead of the structure itself, in general the 



actual structure results in a flatter structure in the quasi-transform domain. 

2.3. Layer segmentation using GTDP 

Once the image is transformed into the quasi-polar domain, it can be segmented using our 
GTDP-based automatic layer segmentation algorithm [6]. In summary, this technique 
represents the image as a graph containing nodes (i.e. image pixels), and edges, which connect 
adjacent pixels. We assign weights to each of the edges based on a priori information about 
the structure to prefer paths through the structure boundaries. Our GTDP technique includes 
automatic endpoint initialization, where an additional column is added to either side of the 
image with minimal vertical edge weights. As we detailed in [6], the start and end nodes can 
then be arbitrarily defined anywhere within these two columns. Finally, any optimization 
method that finds the minimum weighted path from the start node to the end node (e.g. 
Dijkstra's method [43]) can be used to segment the structure. 

2.4. Transformation back into the Cartesian domain 

Once the structure is segmented in the quasi-polar domain, the segmentation information is 
transformed back into the Cartesian domain by applying the inverse of Steps 1 and 3 from 
Section 2.2. 

2.5. Segmentation of subsequent structures 

After the first structure is segmented. Steps 1-3 in Section 2.2 are repeated for each 
subsequent structure (e.g. for each individual cell or cyst). We can achieve this either by 
sequential segmentation of individual features, or by parallelizing the segmentation to make 
this process more efficient. 

3. Implementation for RPE cell segmentation 

AMD is the leading cause of irreversible blindness in Americans older than 60 years [44]. In a 
recent study. Ding et al. showed that anti-amyloid therapy protects against RPE damage and 
vision loss in dinAP0E4 mouse model of AMD [12]. In [12], the average size and density of 
the cells were used as biomarkers for disease progression. To attain these quantitative metrics, 
22,495 cells were analyzed using customized semi-automatic segmentation software. To 
achieve a drastically faster segmentation rate, this section describes an implementation of the 
generalized closed-contour segmentation algorithm from Section 2 to fully automatically 
segment RPE cells in confocal microscopy images of flat-mounted retina of normal and AMD 
mouse models. 



transform the grayscale image into its quasi-polar equivalent, . 



structure is not perfectly flat in . This also implies that a pilot estimate shape closer to the 
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3.1. Data set 

In the following, we describe our fully automatic segmentation method as tailored to the data 
set from Ding et al.'s report published in PNAS [12]. However, the adaption of our method to 
any other data set is straightforward. The data set in [12] included 23 confocal microscopy 
fluorescence images {(/f^), }f^^ e 9^^'''^^''' (0.64 x 0.64 mm^) of central RPE flat mounts 

from 17 mice with approximately 1000 cells per image. RPE cells were stained with a rabbit 
antibody against ZO-1 (40-2200, 1:100; Invitrogen) and imaged on a Nikon Eclipse CI 
microscope. 

3.2. Pilot estimation of cell morphology 

We first normalized the intensities of the image (^if) from [12] s.t. 0<lf <1. We then 

smoothed if using an 1 1 x 1 1 pixel Gaussian filter with standard deviation a = 25 pixels. 
Next, we computed the extended-minima transform [45] that found all minima and suppressed 
those with a depth less than 0.05. To generate the pilot estimate Bf for all cells, we set 

Bf (z) = 1 for all points z lying within the convex hull [46] of each minima. To localize the 
desired cell and its corresponding pilot estimate, we first cropped if and Bf to 201 x 201 
pixel image blocks corresponding to the largest conceivable cell size in our model. Figure 3(a) 
shows the cropped cell image, and Fig. 3(b) shows B^, the cropped pilot estimate image. 

Both images were centered at the z^^f computed for the desired pilot estimate. Since B^ 
contained other cell estimates, we then set all pixels outside the desired estimate in B^ to 
zero. 

3.3. Image transformation into the quasi-polar domain 

We next transformed I^ -^I^^I^ (Fig. 3(a, d, f)) by following Steps 1-3 in Section 2.2 
with the following implementation details: 

■ After generating B^ in Step 1, we found the logical OR combination of b^j^ = b^^ I b^^ 

and set b^^ = b^^ = b^^ , since = 6^^ = 0° = 360°. 

■ For Step 2, we first segmented the boundary in B^ using the GTDP layer segmentation 

technique. We then smoothed the segmentation using a moving average filter with a 
span of 1%. This smoothed cut was set as the function r = f{0) (Fig. 3(c)), which 
provided the shape information necessary to flatten the image. 

■ To generate the edge weights in Section 3.4, we created a threshold image T , where 

5,(z), if5,(z)>mean(5J ^ 
0, otherwise 

and 5^ = wiener 2(7 J (MATLAB notation; The MathWorks, Natick, MA) was 
denoised using the Wiener filter. We then set all connected components in r. 
smaller than 20 pixels to zero, and for all z where Bf {z) = l and B^ (z) = 0, we set 
T^{z) = -l. Finally, we transformed J. into its quasi-polar equivalent, T^. 

While Fig. 3(f) shows the resulting quasi-polar transformed cell image 7^, Fig. 3(e) shows 

the quasi-polar transformed pilot estimate. As can be seen, the transform yielded a fully 
flattened pilot estimate and an approximately flattened cell structure. 
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Fig. 3. RPE cell segmentation using the quasi-polar transform, (a, b) Image containing the cell 
to segment (a), and the pilot estimate (b) with its centroid (blue asterisk), (c) Polar- 
transformation of (b) segmented using GTDP to generate r - f{6) (green), (d) Polar- 
transformation of (a) using the centroid from (b) as the reference pixel. The black regions are 
invalid points that lie outside the image in the Cartesian domain, (e, f) Lnages (c) and (d) 
flattened based on r = f{0), respectively, (g) Transformation of r' = f(0') from (f) back 
into the Cartesian domain. 



3.4. Cell segmentation using GTDP 

Once the quasi-polar transformed image was generated, we used our GTDP layer 

segmentation framework to segment the cell boundary. To do so, we defined the following 
three cost functions: 

(1) Penalize bright nodes further from the centroid to avoid segmenting multiple cells at 

once(ri inEq. (7)). 

(2) Penalize nodes that include the pilot estimates of other cells to also avoid segmenting 

multiple cells at once (T 2 Eq. (7)). 

(3) Penalize nodes that fall below the threshold criteria for T (T 3 in Eq. (7)). 

We implemented these three cost functions using the following MATLAB (The 
MathWorks, Natick, MA) notation, where and t^^ are corresponding columns in T and 

T , respectively, such that 



t^^ = bwlabel(r^ > 0), 



■- bwlabeK-- 1 



-1), 



(6) 



t^^=2' bwlabel(r^ -1) - 2. 
We then combined these cost functions such that 



^.3(^.3<0) = max(r^3) + 2, 



(7) 



^q ^ql^^qS- 



Finally, we calculated the edge weights using 

w^^ = normalize (-/^ (z, ) (z^ ), 0, 0.25) 
+ normalize [t^ (z, ) + (z, ), 1, 2) + 



(8) 
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where w^^ is the weight assigned to the edge connecting pixels and z^, and 
^min =0.00001. The nomialize(x, z) notation indicates a normalization of the values x to 
range from yioz 

Prior to finding the shortest path, we discarded all edges from the graph containing nodes 
in previously segmented cells to avoid segmenting them again. We also excluded the edges 
that allowed those cells to be included within a newly segmented cell. We then found the 
minimum- weighted path p using Dijkstra's method, where r' = p{0') and 0' represent the 

row and column axes of , respectively. Since 0\=0\ in the quasi-polar domain, r \ 
must equal r\\ in practice however, this is not guaranteed when using automatic endpoint 
initialization (Section 2.3). Therefore, if r\^r\, we set both endpoints equal to {r\,6'^ 
and found the corresponding minimum- weighted path p^. We then set both endpoints equal 
to {r\,6\) and found the path p^. The ultimate selected path p (Fig. 3(f), magenta) was 
the one that passed through the brighter overall path, where 

„ |a. ifi:/,{A(^'),^') < E/,(Pe(^',X^',) 
Pq , otherwise 

3.5. Cell transformation back into the Cartesian domain 

We transformed the segmentation r' = p{9') back into the Cartesian domain (Fig. 3(g)) by 

following the method in Section 2.4. However, since the cell was segmented using only a 
cropped section of the image, we shifted the coordinates of the segmented layer back to its 
appropriate location on the original- sized image ^ 9^1024x1024 ^j^^ result was a binary edge 
image g 9^10^4x1024 ^ ^y^q^q pixels {z^ coincided with the segmentation such that 
E^{z^)-l. For all pixels z^ and z^+j that were unconnected, we used linear interpolation to 
connect the dots and create a closed-contour. 

3.6. Iteration for all cells 

To segment all other cells, we iterated the steps outlined in Sections 3.3 through 3.6 for all 
pilot estimates of the cells in the image. To prevent gaps from occurring between adjacent 
cells, we created a preference for already- segmented cell borders by brightening their intensity 

values, where /^(z) = max(/^(z),mean(/^)-Fl.5-std(/^)) and r(z) = l for all z where 
E,{z) = l. 

3. 7. Refinement 

The gold standard segmentation from [12] used in our validation study (Section 4.1) separated 
cells by drawing a line through the middle of the cell border with a single pixel thickness. 
Since the presented algorithm did not necessarily yield a single-pixel cell border, this section 
describes the post-processing steps required to match our automatic segmentation with the 
gold standard protocol, and to remove any erroneous cell segmentations. 

First, we removed all cells smaller than 50 pixels by setting E^{z) = l for all z lying 

within the cell. We then generated the skeleton [47] of E^ and removed all spurs [47]. Next, 

we looked at individual cell edges by setting the branch points in E^ equal to zero. We 

defined the minimum edge threshold to be ^ = mean (/^(Z^))- 1.5- std(/^(Z^)), where 

was the set of all pixels where ^^(z) = 1. If all pixels for a given cell edge fell below t, or if 
80% of the pixels fell below t and there were at least 15 such pixels, then the cell edge was 
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removed from . Once we removed all invalid cell edges, we restored all branch points on 
and removed any newly created spurs or pixels no longer connected to the cell edge. 

4. Experimental results 

This section shows the closed-contour segmentation results obtained using the procedures 
described in Section 3. Section 4.1 evaluates in detail the RPE cell segmentation algorithm 
described in Section 3, while Section 4.2 briefly illustrates the segmentation of closed-contour 
features for other ophthalmic applications. 

4.1. RPE cell segmentation results 

To validate our fully automatic segmentation results against the gold standard, we compared 
the cell count and mean cell size for each of the 23 images in [12]. To prevent from biasing 
our results, these images were not used during algorithm development stage. Instead, a single 
confocal microscopy fluorescence image of an AMD mouse model not included in the 
validation data set was used to train the algorithm. 

We did not alter the data in Ding's study in any shape or form, with one slight exception. 
In [12], closed-counter features smaller than 100 pixels were absorbed into the neighboring 
cells, since such regions were likely a result of incorrect manual segmentation. In our 
replication of the results from [12], we instead discarded these regions from our quantitative 
comparison since such regions were negligible. In [12], regions of the image considered 
invalid due to imaging errors, flat-mount preparation errors, or cells being cut off by image 
borders, were marked separately and discarded from the analysis. In this study, we used the 
markings from [12] to ignore the exact same invalid regions in our analysis. Note that, the 
above two types of outliers occupied a very small percentage of the total area segmented. 

Quantitative results for the RPE cell segmentation algorithm described in Section 3 are 
shown in Table 1 . The average error between automatic segmentation and the gold standard 
was 1.49 ± 1.44% for the cell count and 1.32 + 1.53% for the mean cell area. The median 
error was found to be 0.97% and 0.71% for the detected number of cells and mean cell area, 
respectively. Qualitative results for this validation study are shown in Fig. 4, where Fig. 4(a) 
(corresponding to Image 16 in Table 1) is the automatically segmented confocal fluorescence 
image of an AMD flat-mounted mouse retina, and Figs. 4(b-g) are zoomed-in cell 
segmentation results. Figure 5 shows two confocal images of RPE cells (Figs. 5 (a, d)) along 
with their gold standard (Figs. 5(b, e)) and automatic (Figs. 5(c, f)) segmentation results. The 
image shown in Figs. 5(a-c) corresponds to Image 9 in Table 1 with approximately median 
error, and Figs. 5(d-f) show Image 10 from Table 1 with maximum error. Lastly, Fig. 6 shows 
zoomed-in images of Images 6 and 9 from Table 1 comparing the gold standard (Figs. 6(a-e)) 
to automatic segmentation (Figs. 6(f-j)). The images in Figs. 6(a-c) are examples of erroneous 
gold standard segmentation, while Fig. 6(j) is an example of erroneous automatic 
segmentation. Their corresponding images (Figs. 6(e-h)) are examples of accurate 
segmentation. Finally, Figs. 6(d, i) show an undetermined case where it is unclear which 
segmentation is correct. 

All 23 images used in the study and their corresponding manual and automatic 
segmentation data are available at http://www.duke.edu/~sf59/Chiu_BOE_2012_dataset.htm. 
The "single click" algorithm was coded in MATLAB (The MathWorks, Natick, MA) and 
resulted in an average computation time of 2.95 minutes per image (with an average of 964 
cells per image) on a desktop computer with a 64-bit operating system, Core 17 2600K CPU 
(Intel, Mountain View, CA), solid state hard drive, and 16 GB of RAM. This time included 
the overhead required for reading and writing operations. 

4.2. Segmentation results for other applications 

Section 2 presents a generalized segmentation framework for closed-contour structures, and 
Sections 3 and 4.1 discuss the implementation and results for RPE cell segmentation. We have 
additionally implemented the algorithm to segment other ophthalmic features. These include 
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Table 1. Comparison of the RPE cell count and average cell area obtained for each 
confocal fluorescence microscopy image via fully automatic segmentation versus the gold 

standard 







Number of Cells 






Mean Cell Area 






Gold 






Gold 








Standard 


Automatic 


Error 


Standard 


Automatic 


Error 


Image 


(n) 


(n) 


(%) 


(pixels) 


(pixels) 


(%) 


1 


885 


855 


3.39 


1003 


1036 


3.27 


2 


710 


721 


1.55 


1228 


1204 


1.97 


3 


829 


830 


0.12 


1087 


1082 


0.47 


4 


776 


775 


0.13 


1126 


1124 


0.23 


5 


825 


817 


0.97 


1092 


1100 


0.71 


6 


923 


902 


2.28 


866 


882 


1.83 


7 


981 


937 


4.49 


899 


940 


4.53 


8 


971 


960 


1.13 


875 


881 


0.72 


9 


1097 


1088 


0.82 


832 


836 


0.52 


10 


1253 


1181 


5.75 


715 


757 


5.94 


11 


1187 


1170 


1.43 


771 


781 


1.24 


12 


833 


828 


0.60 


1062 


1065 


0.21 


13 


900 


895 


0.56 


1007 


1010 


0.23 


14 


1235 


1220 


1.21 


730 


736 


0.89 


15 


1005 


999 


0.60 


892 


895 


0.32 


16 


1109 


1075 


3.07 


815 


839 


2.96 


17 


1084 


1077 


0.65 


836 


840 


0.37 


18 


1003 


982 


2.09 


916 


934 


1.93 


19 


1013 


1008 


0.49 


865 


866 


0.11 


20 


1054 


1042 


1.14 


856 


863 


0.90 


21 


931 


927 


0.43 


970 


972 


0.16 


22 


973 


967 


0.62 


931 


935 


0.36 


23 


919 


912 


0.76 


964 


969 


0.45 


Median 


973 


960 


0.97 


899 


934 


0.71 


Mean 


978 


964 


1.49 


928 


937 


1.32 


Std 


141 


132 


1.44 


130 


123 


1.53 





(b) 



(c) 




Fig. 4. Examples of RPE cell segmentation, (a) Automatically segmented confocal 
fluorescence image of a flat-mounted mouse retina (Image 16 in Table 1). The cell borders 
could still be segmented in cases when (b) the pilot estimate (yellow) was off-center and not a 
close estimate of the cell, (c) image artifacts or extraneous features were present, (d) the 
reflectivity of the cell border varied, (e) the cell had a low signal-to-noise ratio, (f) the cell was 
of abnormal shape, and (g) cell sizes were large or small. Each colored box in (a) corresponds 
to the zoomed-in image with the same colored box in (b-g). 
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Gold Standard Automatic Segmentation 



Image 10 




(d) (e) (f) 



Fig. 5. Comparison of fully automatic segmentation versus the gold standard, (a, d) Confocal 
fluorescence microscopy images of flat-mounted mouse retina, (b, e) Gold standard 
segmentation of RPE cells (magenta) obtained semi-automatically using an independent 
technique, (c, f) Fully automatic segmentation (magenta) using the presented closed-contour 
GTDP technique. For the gold standard, cells bordering the image and invalid regions due to 
folding of tissue during preparation and imaging artifacts were ignored. These regions were 
therefore discarded (a-f , black borders) for the comparison study shown in Table 1 . 



Image 9 Image 9 Image 



Gold Standard 





/T 









Automatic 
Segmentation 




Fig. 6. Zoomed-in comparison of fully automatic segmentation versus the gold standard, (a) 
Erroneous gold standard segmentation : The small middle cell was merged with adjacent cells. 

(b) Erroneous gold standard segmentation that did not significantly alter quantitative 
comparison : Although the middle cell was incorrectly shifted, the cell count remained correct. 

(c) Erroneous gold standard segmentation : An enlarged diseased cell was incorrectly separated 
into two cells. We emphasize that such errors (a-c) were very infrequent in the gold standard 
data set consisting of thousands of semi-automatically segmented cells. However, existence of 
even a handful of such errors shows the limitation of subjective segmentation techniques 
relative to the automatic segmentation (f-h). (d, i) An undetermined case : The gold standard (d) 
delineated two separate cells, while the automatic segmentation (i) grouped them as a single 
cell. Since these are diseased RPE cells, it is unclear whether cells with a partially disappeared 
cell border should be considered individually or as a unit, (j) Erroneous automatic 
segmentation : Borders were incorrectly drawn through the cells due to an artifact, (e-h) 
Accurate gold standard (e) and automatic (f-h) segmentation. 
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(a) (b) 

Fig. 7. (a) SDOCT retinal image of a patient with DME. (b) Fully automatically segmented 
retinal layers (inner limiting membrane - blue; RPE - magenta; Bruch's membrane - cyan) via 
layer GTDP, and segmented intra-retinal cysts (magenta) via closed-contour GTDP. 




Fig. 8. (a) AOSO image of cone photoreceptors provided by the authors of [18]. (b) Fully 
automatically segmented cones using the closed-contour GTDP segmentation technique. 



intra-retinal cysts seen on SDOCT images (Spectralis; Heidelberg Engineering, Heidelberg, 
Germany) of patients with DME (Fig. 7). We have also included automatic segmentation of 
cone photoreceptors on AOSO images [18] of the retina (Fig. 8). To segment the cysts and 
photoreceptors, the same basic framework as described in Section 3 was utilized. Changes 
required to adapt the algorithm to different image types included a modification of image 
filters and pilot estimation techniques, changing the edge weights, and different refinement 
steps to exclude erroneous segmentation (e.g. removing segmentations of hypo -reflective 
blood vessels in SDOCT images). 

5. Discussion 

To the best of our knowledge, no other technique has been reported for fully automatically 
segmenting confocal fluorescence images of RPE cells. The qualitative results shown in Fig. 4 
demonstrate accurate automatic cell segmentation despite the presence of image artifacts, low 
signal-to-noise ratio, or inaccurate estimations. Furthermore, the cell shape and size were not 
constraining factors for the presented automatic algorithm. 

Quantitative results show that our fully automatic algorithm accurately segmented RPE 
cells on confocal images of flat-mounted mouse retina with AMD, yielding cell count and 
mean cell area measurements with an average error of less than 1.5%. Automatic 
segmentation errors similar to Fig. 6(j) occurred due to the presence of bright imaging 
artifacts, which erroneously altered values in the first cost function of Eq. (6). However, as 
shown in Figs. 6(a-c), even the gold standard segmentation reported in PNAS contained errors 
as well. Such errors will naturally occur in most manual or semi-automatic segmentation tasks 
required for large-scale studies, where multiple human experts subjectively grade images at 
different dates. Furthermore, since manual correction is extremely time consuming, cells with 
inconsequential errors (Fig. 6(b)) may not have been a priority to correct. As a result, a 0% 
error from the gold standard does not imply perfect segmentation. 

To further reduce the effect of pilot estimation and achieve more accurate segmentation, 
one may employ an iterative approach, in which the result of segmentation for one iteration is 
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used as the pilot estimate for the next iteration. In addition, application of modern adaptive 
denoising methods to raw images, which remove imaging artifacts without altering the 
anatomic structures, may further improve the automated segmentation performance [48,49]. 

Our motivation for using cell count and area as validation metrics stems from the need to 
quantify the extent of cell damage present in different groups of diseased mice [12]. While 
error in the absolute position of each automatically segmented cell may have been a more 
direct validation metric, the presence of errors in the gold standard made this unfeasible. 
Furthermore, the cell boundary is several pixels thick, making it difficult to assign a "true" 
boundary position. While the gold standard divided cells often through the middle of the cell 
border, this was not always the case. As a result, we deemed cell count and area to be a more 
meaningful and viable validation metric. 

Finally, we showed preliminary results demonstrating the segmentation of cysts seen in 
SDOCT images of eyes with DME (Fig. 7) and cone photoreceptors seen in AOSO images of 
the retina (Fig. 8). Prior studies on retinal fluid segmentation have focused on the detection of 
larger regions of intra-retinal fluid [19] or fluid constrained by retinal layer segmentation [21]. 
However, our algorithm is sensitive to the small changes in intra-retinal fluid buildup. With 
regards to cone photoreceptor segmentation, many previous works assumed a given 
photoreceptor size, shape, or packing [32-35,37]. Our method relaxes these constraints, 
allowing for a future version of the algorithm capable of segmenting both photoreceptor rods 
and cones in a single image. Lastly, while we plan to validate these algorithms in future 
studies, we believe these preliminary results are encouraging, as they show the flexibility of 
our algorithm to be extended to different structures, diseases, and imaging modalities. 

6. Conclusion 

In summary, we demonstrated the extension of our automatic GTDP framework to segment 
not only layered, but also closed-contour structures. Implementation of the algorithm for RPE 
cell segmentation in confocal fluorescence images of flat-mounted AMD mouse retina 
resulted in an accurate extraction of cell count and average cell size. We believe that such a 
tool will be extremely useful for future studies, which use cell morphology as a biomarker for 
the onset and progression of disease. While validation has yet to be shown for other 
ophthalmic applications, we have demonstrated preliminary results for intra-retinal cyst 
segmentation in a SDOCT image as well as cone photoreceptors segmentation in an AOSO 
image. This is highly encouraging for reducing the time and manpower required for 
segmenting closed-contour features in ophthalmic studies. 
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